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ABSTRACT 

The formation of the giant planets in our solar system, and likely a majority of giant exoplan- 
ets, is most commonly explained by the accretion of nebular hydrogen and helium onto a large 
core of terrestrial-like composition. The fate of this core has important consequences for the evo- 
lution of the interior structure of the planet. It has recently been shown that H2O, MgO and Si02 
dissolve in liquid metallic hydrogen at high temperature and pressure. In this study, we perform 
ab initio calculations to study the solubility of an innermost metallic core. We find dissolution of 
iron to be strongly favored above 2000 K over the entire pressure range (0.4-4 TPa) considered. 
We compare with and summarize the results for solubilities on other probable core constituents. 
The calculations imply that giant planet cores are in thermodynamic disequilibrium with sur- 
rounding layers, promoting erosion and redistribution of heavy elements. Differences in solubility 
behavior between iron and rock may influence evolution of interiors, particularly for Saturn-mass 
planets. Understanding the distribution of iron and other heavy elements in gas giants may be 
relevant in understanding mass-radius relationships, as well as deviations in transport properties 
from pure hydrogen-helium mixtures. 

Subject headings: planets and satellites: interiors, planets and satellites: dynamical evolution and stabil- 
ity, planets and satellites: individual (Jupiter, Saturn) 



Despite recent advances in computational 
methods improving understandin g of the hydrogen- 
helium dominated outer layers ( McMahon et al 



2012tlFrench et al.ll2012l:lMilitzeij|2013tlWilson k Militzer 



20101 ). knowledge of the deep interior structure of 
giant planets is limited. Determining the size of 
a dense central cores in a giant planet is depen- 
dent upon the model and equation of state used. 
Current observational evidence yields rec ent esti- 
mate s for present day core mass of ~0— 10 (IGuillot 



20081) Earth 



20051 ) and ~14-18 (Militzer et al 
masses for Jupiter, and ^9—22 (jGuillotl 120051 ) 
Earth masses for Saturn. The Juno spacecraft, en 
route to Jupiter, will improve this constraint with 
more precise measurements of the giant's gravita- 



tional field (|Helledll201l[ ). Meanwhile, the density 
profiles of Neptune and Uranus allow non-unique 
solutions for th e compositional struc ture for much 

of the interior (lGuillotlll999bll2005l). 

It has long been suggested ( Stevensonlll982al lb). 
that a portion of this dense material might be 
redistributed in solution with hydrogen. As a 
result, erosion of a dense core would cause it 
to shrink over the lifetime of the planet. Pos- 
sible consequences of this process are only be- 
gin ning to be enumerated in evolutionary mod- 
els (IChabrier fc Baraffd lioir? : iLeconte fc Chabrien 
2012tlMirouhll2012l) . The establishment of a gra- 
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dient in concentration of a heavy dissolved com- 
ponent may change the nature of convection in 
a portion of the planets interior. This 'double- 
diffusive' is hypothesized to reduce the efficiency 
of heat transfer, thereby altering the thermal evo- 
lution of the planet. Comprehensive understand- 
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ing of the process has been hmited by the laek 
of knowledge of the solubihty of various phases in 
metallic hydrogen, as well as poor understanding 
of the sealing of eonveetive efheiency in the pres- 
ence of competing gradients of composition and 
temperature. In this study, we address the first 
issue for iron metal. 

As a result of con tinuing discoveries by Kepler 
(|Borueki et al.ll2010t ) and other exoplanet surveys, 
the number of confirmed planets has climbed to 
over 800, the majority of which are giants. This 
presents a growing sampling of planetary mass- 
radii relationships that will be fundamental to un- 
derstanding the evolution of giant planet interiors. 
The range of mass-radius relationships observed 
for exoplanets exhibit variation beyond those in 
the sola r system. In som e eases, such as Corot- 
20b (Deleuil et al. l l201lh . the relationships may 
even defy explanation by simple structural models. 
Redistribution of dense core material lowers the 
heavy element content required to explain anoma- 
lously high observed densities. 
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Fig. 1. — Gibbs free enerby of solvation for solid 
Fe in liquid metallic hydrogen. Negative values 
favor dissolution for a solute ratio of 1:256. 

The favored model for gas giant forrnation 
(iMizuno et al.ll978HBodenheimer fc PoUackllQSl 



Pollack et al .1 119961) relies on the early formation of 
a large planetary embryo of critical mass to cause 
runaway accretion of hydrogen and helium gas. A 
competing theory involves collap se of a region of 
the disk under self-gravity, e.g. ( Bossl[l997' ). but 
may have difficulty explaini ng significant enrich- 
mcnt of refractory elements ([Hubbard et al ] |2002t 
[Cuillot 2005). The immediate result of a core- 
accretion hypothesis is a planet with the ice-rock- 



metal embryo residing at the center as a dense 
core, surrounded by an extensive layer of metallic 
hydrogen and helium. The role of core erosion to 
the subsequent evolution is a major source of un- 
certainty, but in principle, can explain shrinking 
of cores to masses smaller than those necessary to 
form the planet under the core-accretion hypoth- 
esis. 

Core erosion in giant planets can be addressed 
by determining the solubility of analogous phases. 
Previous studies have considered an icy layer 



of flui d and superionie II2 O ( Wilson fc Militzer 



2012a 

consisting of 
and Si02 



Wilson et all l2013l). and a rocky layer 



MgO 



(Wilson fc Militzeil l2012b[ ) 
^13^), which have 



([Gonzalez et al. 
been shown to separate at relevant conditions 
([Umemotoll2006l ). Assuming the same gross distri- 
bution, elements as terrestrial bodies, the inner- 
most core component would be a dense, metallic 
alloy composed primarily of iron. 

Ah initio random structure searches ( Pickard fc Need^ 
20091) demonstrate that iron remains in a hexag- 
onal close packed (hep) structure remains sta- 
ble up to pressures approaching Jupiter's center, 
^2.3 TPa, at which point it undergoes a phase 
transition to fa ce centered cubic (fee) structure. 
( Stixrud j 2012) demonstrated a gradual decrease 
in this transition pressure wit h temperature . Sim- 
ulations of liquid hydrogen (iMilitzer et al.ll2008 : 
Militzeil [201a iMcMahon et all I2OI2I) undergo a 
gradual transition from molecular to metallic, 
which is complete by ~ 0.4 TPa at low tempera- 
tures. 

Compression of materials with modern exper- 
imental techniques can reach megabar pressures, 
however the pressure-temperature conditions near 



Table 1: Gibbs free energy of solvation for Fe in 



liquid H 


P 


T 


Fe Phase 


AG 


(GPa) 


(K) 




(eV) 


400 


2000 


hep 


-2.2 ± 0.14 


1000 


2000 


hep 


-2.5 ± 0.12 


1000 


2000 


fee 


-2.3 ± 0.13 


4000 


2000 


fee 


-1.71 ± 0.056 


4000 


15000 


fee 


-8.42 ± 0.066 


4000 


20000 


fee 


-11.34 ± 0.078 
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Jupiter's core (>4 TPa) remain inaccessible. Tem- 
peratures in Shockwave experiments climb rapidly 
at high pressures in comparison with planetary 
adiabats, while diamond anvil cell experiments are 
best suited for low temperatures. As a result, sim- 
ulations based on ab initio theories are best suited 
for directly probing the conditions of gas giant in- 
teriors. 

We performed density functional theory molec- 
ular dynamics (DFT-MD) simulations to deter- 
mine the energetics of a dissolution reaction, in 
which solid iron dissolves in pure liquid hydrogen. 
We calculate a Gibbs free energy of solvation: 



AGsoi (Fe : 256H) 



G (H256Fe) 



(1) 



G(H256) + ^G(Fe32) 



where G(H256) and G(Fe32) are the Gibbs free 
energies of a pure hydrogen liquid and solid or 
liquid iron. G (H256Fe) is the Gibbs free energy 
of 1:256 liquid solution of iron in hydrogen. We 
assume that analysis of a single low-concentration 
solution is sufficient to determine the onset of core 
erosion, since the reservoir of metallic hydrogen 
would be much larger than the core. This does not 
rule out non-ideal effects of higher concentrations 
that might exist in a narrow, poorly convecting 
layer at the top of a core. 

Free energy calculations require the determina- 
tion of a contribution from entropy, which is not 
determined from the standard DFT-MD formal- 
ism. To achieve this, we adopt a two step thermo- 
dynamic integration method, as used in previous 



studies (iMorales et all l2009t 



.. _ IWilson fc Militzer 

1201(1 12ni2al lbir The method requires integration 
of the change in Helmholtz free energy over an 
unphysical, yet thermodynamically permissible, 
transformation between two systems governed by 
potentials Ua (ri) and Ub (ri). We define a hybrid 
potential Ux — {1 — X)Ua + XUb, where A is the 
fraction of the potential t/f, (ri). The difference is 
Helmholtz free energy is then given by 



(2a) 



Jo 



dX {Ub (r;) - Ua (ri))A (2b) 



where the averaging is over configurations gener- 
ated during simulations of the system governed by 
the hybrid potential. 



To increase the efficiency, we calculated the 
Helmholtz free energy in two steps each involv- 
ing an integral of the form of Eq. 2. We first 
find AFdft^cI, between the systems governed by 
DFT and classical pair potent ials, which are fi t via 
a force-matching method (Izv ekov et al. l liool . In 
the second step Ai^ci->an, between the pair po- 
tentials and a reference system with an analytic 
solution Fan. The Helmholtz energy for the DFT 
systems, Fdft = Fa„-hAFDFT^ci + AFci^an, may 
then be compared directly. The first step requires 
DFT-MD simulations for a small number of A val- 
ues, while the second uses applies a much faster 
classical Monte Carlo approach at numerous val- 
ues of A to ensure a smooth integration. 
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Fig. 2. — Breakdown of AGsoi into contributions 
from: internal energy, AUsoi, pressure effects, 
PAVsoi, and entropic effects, —TASsoi- Plots 
show variation with (a) temperature at P=4 TPa, 
and (b) pressure at T=2000 K. 

Finding a suitable reference system is essential 
to the method, as it allows AFci_j.dft to be found 
with a small number of number of integration 
steps, and prevents solids from melting or trans- 
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forming to a new structure. For liquid systems, we 
integrate to a non-interacting ideal gas, using clas- 
sical two-body pair potentials as the intermediate 
step. For solid iron, we integrate to an Einstein 
solid, a system of non-interacting, 3-d harmonic 
oscillators. For the intermediate classical solid sys- 
tem we use a 50-50 'mixture' of two-body and one- 
body harmonic oscillator potentials. Spring con- 
stants for the Einstein terms are found from the 
mean-squared displacement of atoms from their 
ideal lattice sites during a DFT-MD simulation. 

All simulations presented here were performed 
using the Vienna ab initio simulation package 
(VASP) (jKressel ll99(J) . VASP uses the DFT 
formalism utilizing pseudopot entials of the pro- 
jector augmented wave type ( Blochl 19941) and 
the exchange-correlation functional of P erdew, 
Burke and Ernzerhof f|Perdew et al.|[l996l ). The 
iron pseudopotential treats a [Mg]3pd^4s^ elec- 
tron configuration as valence states, and a 2x2x2 
grid of k-points is used for all simulations. Sim- 
ulations on hydrogen and the solution were per- 
formed with a 900 eV cutoff energy for the plane 
wave expansion, while a 300 eV cutoff was used 
for iron. A time step of 0.2 fs was used for all liq- 
uid simulations, a 0.5—1.0 fs time step was used 
for high and low temperature iron simulations re- 
spectively. The AGsoi results were confirmed to 
be well-converged with respect to the energy cut- 
off and time step. Prohibitively long simulation 
times required that convergence with respect to 
finer k-point meshes be verified over a subset of 
configurations generated by a simulation with a 
2x2x2 grid. 

Iron simulations assume an hep or fee struc- 
ture within their respective stabilit y regimes 
(IPickard fc Needd 120091: [Stixrudell2012l ). We con- 
firmed Fe to be sohd up to 20000 K at 4 TPa, and 
to be a liquid at temperatures as low as 15000 
K at 1 TPa. We also confirmed that the Gibbs 
free energy favors hep stability over fee at 1 TPa, 
though the difference is negligible for our subse- 
quent analysis of dissolution. We found 32 atom 
supercells to be sufficient for Fe simulation. Finite 
size effects required that we use large 256 atom 
supercells for hydrogen, to which one Fe atom was 
added for the solution. Cubic supercells are used 
for fee and liquid runs. In order to maintain the 
same number of atoms for the hep an orthogonal 
supercell defined the combination of hexagonal 
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Fig. 3. — Shift in AG sol from a system with and 
Fe:H ratio of 1:256 to 1:100 and 1:1000 in the low- 
concentration limit. 



unit cell vectors a,a -I- b,and c. 

Cell volumes at each temperature were de- 
termined by fitting a pressure-volume polytrope 
equation of state to short DFT-MD simulations. 
The resulting DFT pressures were all within 0.1% 
of the target value. Gibbs free energies were 
computed for the three systems, Fe32, H256 and 
H256Fei, using the thermodynamic integration 
method with simulation times of 1.0 ps for H256 
and H256Fei and 2.5—5.0 ps for Fe. II256 and 
H256Fei runs with A = 1 were extended to 4.0 
ns for precise calculations of the internal energy, 
which allows for determination of the entropic 
component of the Helmholtz free energy. The cal- 
culated energies and entropy are presented in Tab. 
I, along with the density. The Gibbs free energy 
of solvation, calculated using Eq. 1, is presented 
in Tab. II for each pressure-temperature condi- 
tion. A negative AGsoi implies that the Gibbs free 
energy of the solution is lower than that of the sep- 
arated phases. Therefore, dissolution is favored at 
a solute concentration higher than 1:256. 

We find dissolution of iron to be strongly fa- 
vorable at conditions corresponding to the in- 
teriors of gas giants. Fig. 1 shows the vari- 
ation of AGsoi with temperature and pressure. 
AGsoi exceeds —10 eV per iron atom for plau- 
sible temperatures of Jupiter's core. Dissolution 
remains favorable even at temperatures far be- 
low those predicted by model a diabats (jMilitzer 
2013t iMiUtzer fc HubbardI |2013| ). The energet- 
ics are only weakly dependent on pressure, and 
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AG sol becomes increasingly negative with decreas- 
ing pressure. The solubiUty increases with a nearly 
linear trend in T, yielding slope of ~0.53 meV/K. 
As a result, solubility is favored through the en- 
tire range of conditions considered, and likely the 
entire range for metallic hydrogen regions of giant 
planets. 

Fig. 2 shows a breakdown of the data into con- 
tributions by various thermodynamic parameters. 
Included in the figure are: AFgoi, AUsoi, PAVsoi 
and —T AS sol, respectively, the Helmholtz free en- 
ergy, internal energy, volumetric work and cn- 
tropic contributions contributions to AG sol- Note 
that AFsoi, AUsoi and PAVsoi are calculated in- 
dependently, while AGsoi = AFsoi + PAVsoi and 
—TASsoi = AFsoi — AUsoi are derived. The trend 
of solubility with temperature is dominated by the 
entropic term. The high solubility at low temper- 
atures is reflected in the negative values of AUsoi, 
indicating that the mixed system is energetically 
favorable independently of the entropy term. 

Our calculations neglect any interactions be- 
tween iron atoms in solution, and thus represent 
solubility in the low-concentration limit. AGsoi 
can be related to the volume change associated 
with the insertion of an iron of atom into hydro- 
gen, as other contributions are constant with re- 
spect concentration. It can be shown that results 
for simulations with a l:n solute ratio can be gen- 
eralized to a ratio of l:m using 



tional numerical challenges c ompared to other 



AGc 
knT 



FoiHraFe) - Fo{Hm) - FoiFe) (3a) 
- [FoiH^Fe) - FoiHrn) - Fo{Fe)] 

j [F(H„Fe)-|-is_»mO]'"^' l 
_ f [V(H„Fe)]"^M 



(3b) 



where AGc ~ AGsoi{^ ■ fn) — AGsoi{^ '■ n), and 
F(II„) and V^(II„Fe) are the volumes for the simu- 
lations of hydrogen and the solution respectively. 
Fig. 3 shows the shift of AGsoi at 4 TPa at Fe 
concentrations of 1:100 and 1:1000. AGsoi is de- 
creased for higher concentrations, but not to an 
extent where dissolution would become disfavored. 
At 20000 K this difference between 1:100 and 1:256 
is <2 eV per iron atom, and at 2000 K is smaller 
than uncertainty in calculated values of AGsoi ■ 
The nature of the Fe-H system poses addi- 



solut es considere d previously ( Wilson &: Militzer 



2010i [2012ab: Gonzalez et al 



20131 ). These can 



largely be attributed to the comparatively large 
change in volume and electron density associated 
with the insertion of an iron atom into metallic 
hydrogen. We found it more efhcient to deter- 
mine cell volumes by fitting an equation of state 
to a collection of MD simulations at constant vol- 
ume, rather than performing extended constant 
pressure simulations. Finite size effects were also 
found to be more significant for the Fe-H system, 
due to iron's relatively large volume and number 
of valence electrons. Fig. 4 shows the conver- 
gence of a difference in internal energy between 
H and H-Fe for MD simulations with 128, 256 
and 512 hydrogen atoms. We find 256 hydrogen 
atoms to be necessary in contrast to the previ- 
ous studies that required only 128. The conver- 
gence with k-point grid resolution is also slower 
than in previous studies, and presents the great- 
est uncertainty in this study. MD calculations 
with a 3x3x3 k-point grid are prohibitively ex- 
pensive. An estimate of this error for the re- 
sults presented here is obtained by evaluating 
the internal energy over configurations sampled 
from a MD trajectory with a 2x2x2 k-point grid. 
Fig. 5 shows this estimated correction of AGsoi 
for the k-point grid used. The shifts are on the 
order '^l eV per iron atom, but are consistently 
negative for both quantities, leading to dissolution 
being more favorable. 
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Fig. 4. — Energy of insertion for a single Fe atom 
into supercells containing 128, 256, and 512 atoms. 
Finite size effects are significant for H128, but are 
negligible within error for H256. 
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Wit h the results of previou s studies (jWilson fc Mihtzer 
2012al lhl: iGonzalez et al we can now present 

a comprehensive picture for the solubihty of all 
typical core materials in liquid metallic hydrogen. 
Dissolution is strongly favored for both iron and 
water ice. However, for water the high solubil- 
ity is attributed entirely to the entropy, whereas 
iron has a favorable internal energy component 
that favors dissolution at low temperatures. Both 
phases are found to be soluble throughout the 
entire metallic hydrogen region of both Jupiter 
and Saturn. The rocky components, MgO and 
Si02, have more moderate solubilities, with Si02 
being slightly higher. The saturation curves are, 
however, less steep than the adiabats for Jupiter 
and Saturn. As a result, solubility is favored for 
Jupiter's core, but the rocky components of Sat- 
urn's core may be stable given the present uncer- 
tainty in the planet's adiabat. 
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Fig. 5. — Estimated corrections to AGsoi and 
AGsoi coarseness of k-point grid used in DFT-MD 
runs. DFT calculations with a 3 x 3 x 3 k-point 
grid were performed sampling a trajectory gener- 
ated by an MD simulation with a 2 x 2 x 2 k-point 
mesh. 

The rocky components are found to be sta- 
ble at lower pressures, approaching the metallic 
transition. If Mg and Si were convected upwards 
with sufficient concentration, they may precipi- 
tate, while Fe and H2O would remain in solution, 
at least to the molecular-metallic transition. The 
presence of a significant dissolved component at 
shallow depths may have consequences for the den- 
sity profile and transport properties of hydrogen, 
which influence thermal structure and magnetic 
field generation. 



Core erosion is thermodynamically favorable in 
gas giant planets, with the possible exception of 
smaller, cooler planets, like Saturn. For these 
planets, the outer icy layers are soluble, but the 
rocky layers may not be. The innermost iron com- 
ponent, though soluble, would be isolated from re- 
action with hydrogen. This might allow Saturn 
to have a larger, less eroded core than Jupiter, a 
result consistent with current observational con- 
straints. Nevertheless, our results imply that con- 
firmation of a massive core for Jupiter would sup- 
port the core-accretion model over gravitational 
collapse. While erosion of such a core may be 
slow due to ineflncien t doubl e-diffusive c onvec- 
tion (IStevensonlll982at IChabric r & Baraffe 12007 : 
Leconte fc Chabrieij|2012t lMirouhll2012t) . settling 
of dispersed refractory material to form a core is 
inconsistent with our results. Late formation of a 
core would require a large amount material from 
captured planetesimals surviving descent to the 
planet's center. 

It may be possible to attribute some emerging 
trends in exoplanet mass-radius relationships to 
the difference in solubilities between rock and ice, 
or rock and metal. However, as we have shown, 
such thermodynamic differences are likely to only 
be significant in smaller, cooler planets, where re- 
distribution of dense material by double-diffuse 
convection would be least efficient. The energetics 
of the dissolution reaction should be insignificant 
compared to the role of density in the redistribu- 
tion of dense material. The work required to raise 
an iron of atom to the molecular-metallic transi- 
tion is on the order of 1000 eV, whereas the con- 
tribution from the dissolution reaction is — 10 
eV. We conclude that the process of core erosion is 
thermodynamically consistent with ab initio sim- 
ulations of the relevant materials, and its signifi- 
cance warrants close consideration in future mod- 
els of giant planet evolution. 

This work has been supported by NASA and 
NSF. Computational resources at NAS and NCCS 
were used. We also acknowledge Dave Stevenson 
and Felipe Gonzalez for helpful discussions. 
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Table 2 

Thermodynamic parameters derived from DFT-MD simulations. 



p 


T 


System 


Phase 


P 


F 


U 


G 


S 


(GPa) 


(K) 






(g/cTO?) 


(eV) 


(eV) 


(eV) 


(kfc/K) 


400 


2000 


Fe32 


hep 


14.408 


-190.8(4) 


-154.4(0) 


323.3(4) 


211.(4) 






H256 


liq 


1.2709 


-415.6(5) 


-228.(0) 


419.37(9) 


1088.(2) 








liq 


1.5192 


-423.8(1) 


-234.(0) 


427.1(8) 


1101.(6) 


1000 


2000 


Fe32 


hep 


18.279 


-12.4(0) 


18.1(2) 


1000.8(5) 


177.(1) 






H256 


i!q 


1.8916 


23.0(3) 


175.9(6) 


1425.6(6) 


887.(3) 






H256F0 


liq 


2.2534 


20.4(1) 


175.(2) 


1454.7(1) 


898.(3) 


1000 


2000 


FC32 


fee 


18.269 


-10.35(0) 


19.8(0) 


1003.4(6) 


174.(9) 






H256 


liq 


1.8916 


23.0(3) 


175.9(6) 


1425.6(6) 


887.(3) 






H256Fe 


liq 


2.2534 


20.4(1) 


175.2(3) 


1454.7(1) 


898.(3) 


4000 


2000 


Fe32 


fee 


28.374 


754.91(4) 


777.6(2) 


3365.97(1) 


131.(8) 






H256 


liq 


3.6375 


1392.3(8) 


1487.6(4) 


4310.0(0) 


552.(7) 






H256Fe 


liq 


4.3078 


1412.3(8) 


1508.(4) 


4413.4(7) 


55(7). 


4000 


15000 


Fe32 


fee 


27.826 


382.1(9) 


87(5). 


3045.7(0) 


38(1). 






H256 


liq 


3.3618 


-392.9(3) 


1917.(3) 


2763.9(8) 


1787.(3) 






HsseFe 


liq 


3.9865 


-392.2(5) 


194(3). 


2850.7(4) 


1807.(2) 


4000 


20000 


FC32 


fee 


27.550 


176.(9) 


92(2). 


2866.(8) 


43(2). 






H256 


liq 


3.2731 


-1284.5(8) 


208(1). 


1957.2(6) 


1952.(6) 






H256FC 


liq 


3.8824 


-1294.0(9) 


210(4). 


2035.5(1) 


1972.(0) 



8 



